CHAPITRE 12 


ÉQUATIONS AUX DÉRIVÉES PARTIELLES 


Il existe de nombreuses méthodes numériques de résolution des équations 
aux dérivées partielles (Ralston et Wilf, D. D. McCracken, A. D. Booth, 
J. Girerd, E. Durand, voir bibliographie). 

Nous exposerons ici la méthode des différences finies en l’appliquant à la 
résolution de l’équation de Laplace, mais le lecteur pourra reconduire cette 
méthode pour résoudre d’autres types d’équations telles que l’équation de 
Poisson, l’équation des cordes vibrantes, l'équation de diffusion... 


12.1. — ÉQUATIONS AUX DIFFÉRENCES FINIES 


Ayant à résoudre une équation de deux variables du type : 


DV  D2V 
PS eh 
dx dy” 


on peut quadriller le plan x, y par des droites : x = constante et y = constante, 
chaque nœud du réseau ainsi obtenu étant repéré par deux indices I et J. 


Dar 


Fic. 49. 


En un point de coordonnées x, y, on peut spider le la formule de Taylor 
arrêtée au second ordre, soit : 


dV 


: dV 
VC, }) = Vo; Yo) + 3x dx + 3 “dy 


92 
: uns 
d 2 


1 
À dx 0 


ÉQUATIONS AUX DÉRIVÉES PARTIELLES . 193 


Si on applique cette formule aux quatre points du réseau ci-dessus, il vient : 
OV  hl?:a DV 
RE + Sms 
dXo 2 dx 
OV, Ha OV 
dXo 2 dx 
DV  kl?-a D°V 
ane + e es 
dYo 2 dyo 
OV  K2*:a° DV 
PR + . ss 
dYo 2 dyo 


V(xo + hl:a, yo) = V(xo; Yo) + hl'a: 


V(xo — h2°a, yo) = Vo; Jo) — h2'a: 


| VC; Yo La kl:a) — V(Xo: Yo) + kl'a: 


Vo; Yo — k2*a) = V(xo, Yo) — k2:a: 


On peut multiplier la première équation par : 


la deuxième équation par : 


1 
h2-(h1 + h2) 
la troisième équation par : 
1 
k1-(k1 + k2) 
et la quatrième équation par : 
1 
K2:(k1 + Kk2) 


On obtient ainsi un nouveau système de quatre équations, on peut addi- 
tionner les deux premières entre elles et les deux dernières entre elles puis les 
regrouper obtenant ainsi : 

a? [av 2°V 
2. pv 2 
.2 [oxs dy5 
1e V(xo+hl:a, yo) rs V(xo—h2"a, yo) 
h1 h2 
V(xos Yo+kl'a)  V(xo, Yo —k2*a) Il LP: 
| ki ds 2 — m0 * 4e] VGoY0) 


ou en repassant en notations indicielles : 


PL ON 2,41 pes ere) 

0x$ Oypg a (hl+h2l hi h2 
| 1 VA+1,9) VG-17] [1 1: 
+re* | A pa À es | ed, 
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DV 9? 
L’équation de Laplace —; + —; = 0 prend donc la forme : 


| dx?  Oy 
1 À 
re F Tel Te 
21 [VID , VO, 1. po : ea) 
” h1+h2 h1 h2 k1+k2 k1 . K2 
Si on suppose 41 = h2 = het k1 = k2 = k; il vient : 
h° h? h? 
2-(1 +2) VOD = VE I+D+VET- D + VAL D + 5 VOL) 
et si on suppose encore h = k; il vient : 
4 V(L, D) = VO, J + 1) + VO, J — 1) + VOL + 1,9) + VC — 1,0) 


On voit dans ce calcul que les dérivées impaires disparaissent car il y a 
alternance de + a et de — a. Les équations ci-dessus sont donc équivalentes à 
l’équation de Laplace à des dérivées d’ordre 4 près. 


12.2. — MÉTHODE DE RELAXATION 


Si on écrit les équations aux différences finies pour chaque nœud du réseau, 
on obtient un système d’équations : | 


Pure 1 ]. 2h12 [VGI+1) VAI-1) 
2"hI (5 + ce "EDS Ex | He À D. | 


2-h1-h2 a+ L3) , VG-1, 2] 


k1 + k2 k1 | k2 


Comme les valeurs de V sont connues sur les frontières, le système d’équa- 
tions ci-dessus à un nombre d’inconnues égal au nombre d’équations. C’est 
un système linéaire que l’on peut résoudre en utilisant le sous-programme 
SYLEQ étudié au paragraphe 6.1. Généralement le nombre de nœuds du réseau 
est important et le temps de calcul risque de devenir prohibitif, aussi préfère-t-on 
utiliser une méthode par approximation. 

On effectue une première approximation de la fonction V en chaque nœud 
du réseau. Les valeurs prises au départ sont quelconques, mais on convergèra 
d’autant plus rapidement que l’on choisira des valeurs proches de la solution 
du système. | 

On choisit alors un ordre de parcours du réseau, cet ordre est arbitraire 
mais à chaque nouvelle itération on devra parcourir le réseau dans l’ordre 
choisi initialement. | 
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En chaque nœud (I, D), on calcule un résidu : 


2-h1-h2 park h ÈS) 4 2'h1-h2 pee n V2] 


RES= 55 | hi D) 4e | A 2 


1 1 | 
NE aus] En 


—2.h1-h2: | 
La nouvelle valeur de la fonction V en ce point est recalculée de manière 
à annuler le résidu en ce point, soit : 


oo ia 
VE 3) = VE D + RES [2 Ar Re Gr dé re)| 


On passe alors au nœud suivant dans l’ordre de parcours choisi et on cal- 
cule le résidu en ce point à partir des nouvelles valeurs de la fonction V aux 
nœuds encadrant ce point. | 

Lorsque l’on a parcouru ainsi tout le réseau, on recommence au premier 
nœud du réseau et on calcule de nouveaux résidus et de nouvelles valeurs de V. 
On démontre qu’en chaque nœud du réseau, on obtient ainsi une suite qui est 
convergente et dont la limite est la solution du système. 

S’étant fixé une précision, on s’arrêtera donc lorsque tous les résidus cor- 


_ respondant à une itération seront inférieurs à cette précision. 


Remarquons que si on suppose hi = h2 = h et k1 = k2 = k, le résidu 


_ s'exprime : 


+ 5 2 
- RES=V(I J+1)+V(L J—-1)+ 2 LVG+ 1, P+V({-1,)]-2: (+ n) ve J) 


et la nouvelle valeur de la fonction V : 


V', D) = VE, D) + RES] L2 ( d n)| 


enfin, dans le cas où h = K, on a : 


RES = VOL 3 + 1) + VOL 3 — 1) + VOL + 1,9) + VO — 1,1) — 4 VOL) 


D ct : 


V{L, 3) = VA, J) + RES/4 


12.3. — MÉTHODE DE SURRELAXATION 


On peut accélérer la convergence de la méthode en introduisant dans les 
relations ci-dessus un facteur d’accélération. 


UTILE 
LL + 
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C'est-à-dire que la relation : 


! 1 f 
VE D = VE D + Res [2 RAA Gr PRE TR) 


est remplacée par la relation : 


! | . 
V'D = VOD + © Res [eu h1-h2- Gr h2 ki: 0) 


Young (voir l'c ouvrage de Ralston et Wilf en bibliographie) a montré que 
si, pour une enceinte rectangulaire, I varie de 1 à M + 1 et J de 1 à N + 1,on 
peut prendre pour valeur de w 


œ = 1 + 


12.4, — AMÉLIORATION POSSIBLE DE LA PRÉCISION 


En un point de coordonnées X, y, on peut appliquer la formule de Taylor 
arrêtée au sixième ordre, soit : 


dV 
V=Vo+(x—xo): _ +0 yo): 3 Æ 


2 2 | 2 
| CSS Eee +2*(x—x0) (y —70) se +70) =: = 


2,1 dXo'dYo dy 


He ee de) 4320) 7 


d°V 
2x5 dYo 
e °V 


1 DV. 
a [exo +6 6-x ve 
dx 


d°V 
+ (yo) —| + 
(y — 0) — 


DV 
dx , dYo 


DV 
+70) —| + 
(—70) a | 


21 +6-(x—x0)*"(y—70) 


. TV DV 
a +4 x —%0) (yo): RTE 
dx6" dY6 dXo'dY5 


1 5 d°V 4 d°V ; 3 2 
— + [x] — HG = x0) "= Yo)" +10:(x—x0)" (y 70)" 
[ 0) 2x ( o) CO Yo) TR ( o) *(Y—Yo) 
25V a. DV 
— +5 (x—-x,) (y — CREER ES, 
ve A à CAE vue er 


D°V 
+70)". + 
(© —Y0) al 


d°V 5 à 
 ——— + 10 (X—X . = . 
dx3-0y ( 0) “(Y—Yo) 


MEME NET Nc nu o Ne DE 
7 Ua NUL ART 4 ï 
| 4 h ; 
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1 6 XV : dSV k à 
2, 5 X—X e —— +6: X—X . — oo —— +15. X— X ° = s 
7 [ 0) DS (x—xo) ‘(9 —Y0) DEN @—xo) "(y 70) 


SV , 3. OV a ” 
—— +20: (x—-x ———— +15 (x—x0)" "(y 
5 2 ( 0) (Yo) ACTE ( o) *(Y—Yo) 


>SV s OV « 2SV 
+ +6 (x—Xx0) (70) * —© +(7- Pl: 
Dx2.03S ( 0) (Y—Yo) xs (—Y0 Ds 


On peut l’appliquer aux huit points entourant le point xo, yo : 


q_—_——— 


Au nœud I,J+1, on a : X—xo= et y—yo= 
Au nœud 1+1,J, on a : x—xo=0 et y—yo=a 
Au nœud I,J— 1, on a : X—xo=—a et y—yo= 
Au nœud I—1,J, on a : x—xo=0 et y—yo=—a 


Au nœud 1+1,J+1,ona: x—xo=a et y—yo=a 
Au nœud 1+1,J—1,ona: x—xo——a et y—yo=a 
Au nœud I—1,J—1,ona: x—xo——a et y—yo——a 
Au nœud 1—1,J+1,0ona: x—xo=a et y—yo—=—a 


On peut donc écrire : 


VO, J + 1) + VL + 1, 9) + VO, 3 — 1) + VA — 1,3) = 


D2V D2VT at [D*V DVI 2-4 E dSV 
== 6 L——— es ms + 
4 VE) + à +57) 12 ÉRIEE 2x° | 
et : 
VE+LI+D+VE+ LI = D+VE-LJ-D# VOL, J+1) 
4? [a V DV] 4“ ea DV d*V | 
uen rt] ar Lo op *oroy 
4a$ [o$V  9°V SV dSV 
“| dép Te — 
ds 6! DE s T5 bi op nr 


Si on multiplie la première équation par le nombre 4 et si on l’ajoute à la 
deuxième, il vient : 


4 [VO + D+VO+L DEV 1 D+ V1, DJ+ V(+1,14+1)+V(I+1,7-1) 


d?V | 
+ 


- à 20-V(I, J)+6- 
+EV(I-1,J—-1)+V({-1,7+1)=20" V(I, J)+6°a? É 3y? 


Let tp 
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FE : ui 42, El DV DV 
dx*  dy* 12. 1x 
ee | 4-4 
+ Æ ——* 
dx$  dy° 6! 
2 à D?V 
dx? dy? 
par rapport à x et y, on montre facilement que cette équation se réduit à : 


Partant de la relation : = 0 que l’on peut dériver plusieurs fois : 


4.[VL, J + 1) + VA + 1,7) + VOL J — 1) + VA — 1,9) + VA + 1,3 +1) 
+ VE + 1,3 — 1) + VA — 1,3 — 1) + VA — 1,7 + 1) = 20 VOL). 


ou enfin : 


4.V(L D) = - [VE 341) + VA +1, 3) + VOL 11) + V1, D] 


+ : [V(+1, 341) + VA +1, 1—1)+ V(I—1, J—1)+ V(I—1, J+1)] 


Du fait que les termes impairs s’annulent, les équations ci-dessus sont donc 
équivalentes à l’équation de Laplace à des dérivées d’ordre 8 près. 


12.5, — DÉTERMINATION DU PAS DE QUADRILLAGE 


Il est bien évident que la formule de Taylor arrêtée au second ordre sera 
d’autant plus précise que le pas : a du quadrillage sera faible. Mais, par ailleurs, 
plus ce pas sera faible et plus l’occupation mémoire et le temps calcul seront 
grands. 

Suivant que l’on applique la formule à quatre points ou la formule à huit 
points, le premier terme négligé est une puissance quatrième ou une puissance 
huitième du pas de quadrillage. 

Ayant calculé la solution de l’équation de Laplace avec un quadrillage de 
pas : a en appliquant la formule à quatre points; pour savoir si le pas choisi 
est assez fin, on recalcule la solution avec le même quadrillage mais en utilisant 
la formule à huit points. On regarde alors si les résultats sont différents à la 
précision que l’on s’est fixée, sinon on passe à un réseau plus fin. 


12.6. — APPLICATION A LA DÉTERMINATION DU POTENTIEL 
A L'INTÉRIEUR D'UN COAXIAL RECTANGULAIRE 


Nous allons appliquer les principes précédents à la détermination du DAIERS 
tiel à l’intérieur du coaxial nl suivant (fig. S1) : 
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B/A = 0,5 
CJA = 0,5 
D/B = 0,5 


Le conducteur extérieur étant porté au potentiel + 1 et le conducteur 


- intérieur au potentiel O. 


Le potentiel à l’intérieur du coaxial obéit à l’équation de Laplace, on peut 
donc quadriller l’espace inter-conducteurs comme représenté sur la figure 52 


L . et écrire un programme Fortran de résolution. 


DIMENSIOGN V(13,25) 
READ(IL,101)P,PV,SMEGA 
FORMAT(F16.7) 
READ(IL,102)MAXIT 
FOR MAT(I6) 
WRITE(IW,103) 
FORMAT(// B/A = 0.5’) 
WRITE(IW,104) 
FORMAT( C/A = 0.5’) 
WRITE(IW,105) 
FORMAT( D/B = 0.5///) 
Vi =.0 : 

V2 = 1. 

DS 1 J—1,25 

VOL) = V2. 

V(13,7) = V2 

DS 2 I=—2,12 

V(L1) = V2 

V(,25) = V2 

DS 3 J=—7,19 

V(A,7) = VI 

V(10,7) = VI 

DS 4 I=5,9 

V(L7) = VI 

V(L,19) = VI 

DS 50 J=—2,24 

VG,9) = (V2—V1)/3. 
V(2,7) = 2.(V2—V1)/3. 
V(L1,) = (V2—V1)/3. 


1V(12,3) = 2.*(V2-—V1)/3. 
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Œ 
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« 
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nm 
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3 
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FIG. 52. 
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CONTINUE 

DS 51 1=4,10 

V(L6) = (V2—V1)/6. 

VAS) = 2.4(V2—V1)/6. 

V(L,4) = 3.(V2—V1)/6. 

V(L,3) = 4.*(V2—V1)/6. 

V(L2) = 5.*(V2—V1)/6. 

V(I,20) = (V2—V1)/6. 

V(L2i) = 2.*(V2—V1)/6. 

V(L,22) = 3.(V2—V1)/6. 

V(,23) = 4.*(V2—V1)/6. 

V(L,24) = 5.(V2—V1)/6. 

CONTINUE 

ITN = 1 

RESID = .0 

DO 5 I=—2,12 

DO 5 J=—2,24 

IF(1— 4)6,7,7 

IF(I— 10)8,8,6 

IE(J—7)6,9,9 

IE(J—19)5,5,6 

RES = V(—1,))+V(+1,-+V(LJ—1)+V, QUE A #V(L,J) 

IF(ABS(RES)—RESID)13,13,12 

RESID — ABS(RES) 

V(LD) = V(LI)+GMEGA*RES/4. 

CONTINUE 

IF(RESID — P)60,61,61 

ITN = ITN +1 

IF(ITN— MAXIT)10,10,62 

WRITE(IW,106)MAXIT 

[FSRMAT(/ SN NE CONVERGE PAS APRES/16 RELAXATIONS//) 

GO TO 70 

D 75 I—2,12 

DO 75 J—2,24 

IF(1—4)76,77,77 

IF(1— 10)78,78,76 

IF(J—7)76,79,79 

IF(I— 19)75,75,76 

VAUX = V(LI) 

VA,3)=0.2*(V(LJ+1)+V(L41,9)+V(I—1)+V(L—1,7))+-0.05*(V(I+-1,7+ 1 
2 +VA+1,3—1)+V(—1,7—1)+V(A—1,3+1)) 

IF(ABS(V(LJ)—VAUX)—PV)15,75,151 

CONTINUE 

GO TO 152 

WRITE(IW,153) 

FORMAT(//’ IL FAUT UTILISER UN RESEAU PLUS FIN’/) 

GO TO 70 


Il suffit alors d’écrire, en l'instruction 152, un ordre d’impression 


des 


potentiels et en fin de programme, une instruction PAUSE numé- 


rotée : 70. 
Remarquons que l’on peut réduire l’occupation mémoire et le temps de 
calcul en tenant compte de la symétrie par rapport à la ligne : I = 7. 


PELLETIER. 14 
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On a : 
V(8, 3) = V(6, J) 


et l’équation de Laplace pour la ligne : I — 7 peut s’écrire : 
4 V(,3) = V(,J — D) + V(7,J + D) + 2:V(6,7) 


Donc le résidu pour les points : (7,2), (7,3), 7,4), (7,5), (7, (7, 20), (7,21), (7,22), 
(7,23), (7,24) peut se calculer par : 


RES = V(LJ — 1) + VOL 3 + 1) + 2- VA — 1,J) — 4 VOD) 


12.7, — DÉTERMINATION DES CARACTÉRISTIQUES 
D'UNE LIGNE DE TRANSMISSION 
A PARTIR DE LA DISTRIBUTION DE POTENTIEL 


Connaissant la distribution du potentiel à l’intérieur d’une ligne formée 
de deux conducteurs parallèles à un axe oz et parcourue par des ondes planes, 
on peut déterminer diverses caractéristiques de la ligne par l’application du 
théorème de Gauss (H. E. Green, W. S. Metcalf, voir bibliographie). | 

Ayant quadrillé comme précédemment l’espace inter-conducteurs, on peut 
entourer la section du conducteur intérieur par un rectangle passant par les 
centres de certains carrés du réseau. On peut alors déterminer le flux ® du 
vecteur Déplacement à travers la surface de hauteur unité selon l’axe oz et 
ayant pour base le rectangle défini ci-dessus. 


| 
® 
| 
| 
Û 


Le flux d®,, (fig. 53) à travers la surface rectangulaire de base af et de 
hauteur Az = 1 est : | 
d®,g = — egrad V:dS 


oùeest la permittivité du milieu. 
Soit : 


ui 
4 = | — &° > -dx 


: WE VH) 


EH + 


d®,y = — 


